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Abstract 

A variational formulation of the time-dependent linear response based on 
the Sternheimer method is developed in order to make practical ah initio 
calculations of dynamical spin susceptibilities of solids. Using gradient den- 
sity functional and a muffin-tin-orbital representation, the efficiency of the 
approach is demonstrated by applications to selected magnetic and strongly 
paramagnetic metals. The results are found to be consistent with experiment 
and are compared with previous theoretical calculations. 
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Full wave-vector and frequency dependent spin susceptibility x is a central quantity 
in understanding spin fluctuational spectra of solids. Its knowledge accessible directly via 
neutron-scattering measurements is important due to significant influence of spin fluctu- 
ations to many physical properties and phenomena such, e.g., as the electronic spe- 
cific heat, electrical and thermal resistivity, suppression of superconductivity for singlet 
spin pairing, etc . In magnetically ordered materials, transverse spin fluctuations are spin 
waves whose energies and lifetimes are seen in the the structure of transverse susceptibility. 
High-temperature superconductivity, a highly exciting phenomenon, whose origin is still not 
recognized, can be due to a spin fluctuational mechanism [Q]. 

Despite large past efforts put on the development of methods for ah initio calculations of 
the dynamical spin susceptibility based either on the random-phase-approximation (RPA) 
decoupling of the Bethe-Salpeter equation [Q, or within density functional formalism 
quantitative estimates of x with realistic energy bands, wave functions, and self-consistently 
screened electron-electron matrix elements are scarce in the literature [§-0. This is not only 
due to the absence of complete theory for the proper description of exchange-correlation ef- 
fects which is a true many-body problem, but also because standard perturbative treatment 
of an electronic response has serious problems connected with the summation over high- 
energy states and matrix inversion. 

This paper proposes a method which avoids the latter two problems. The method is 
a time-dependent generalization of an all-electron Sternheimer approach [|l^ which has 
been proved to be very efficient in ab initio calculations of phonon dispersions, electron- 
phonon interactions and transport properties of transition-metal materials including high- 



Tc superconductors |^ . The method employs a muffin-tin-orbital representation [|I2| which 
allows to greatly facilitate the treatment of localized states such, e.g., as d- and f-electrons 
of strongly paramagnetic and magnetic materials whose studying is the main purpose of this 
work. 

Applications to transverse spin fluctuations in Fe and Ni as well as calculations of para- 
magnetic response in Cr and Pd demonstrate an efficiency of the approach and resolve some 
discrepancies found in previous theoretical studies. In particular, experimental evidence of 
an optical spin-wave branch for Ni and its absence for Fe |]T^ is correctly described by 



the present calculation which was not done in either early semiempirical approaches 



or within a recent frozen-magnon scheme |T3]. For the first time, the dynamical suscepti 



bility is calculated ah initio for paramagnetic Cr, a highly interesting material due to its 
incommensurate antiferromagnetism |TB[ . The calculation predicts a wave vector of the spin 



density wave (SDW), and clarifies the role of Fermi-surface nesting. Strong long-wavelength 
spin fluctuations of Pd are evident from the present and earlier theoretical studies. 
The description of the method starts by considering a small external magnetic field 

5^ext{rt) = 5be'('i+^)'"e'"*e-''l*l + c.c (1) 

applied to a solid. Here 6h = J^ii Sb'^^fj. shows a polarization of the field {/i runs over x, y, z 
or over —1,0,1), wave vector q lies in the first Brillouin zone, G is a reciprocal lattice 
vector, and rj is an infinitesimal positive quantity. If the unperturbed system is described 
by charge density p(r) and, in general, by magnetization m(r), the main problem is to find 
self-consistently first-order changes Sp{rt) and (5m(rt) = X^i/ ^^i/(r^)e'^ induced by the field 
5Bea;t(rt). If the polarization 6h in (|l]) is fixed to a particular /ith direction, and 6m{rt) is 
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calculated afterwards, a /itli column of the spin susceptibility matrix ^,0;) will 



be found [17|. This essentially solves the problem. 

A central issue of employing time-dependent (TD) density functional theory (DFT) |ll 
to find the quantities Sp{rt) and 6m{rt) is now discussed. The unperturbed density and 
magnetization are described accurately by the static DFT and are expressed via occupied 
Kohn-Sham states. This is by now a well established method in practical ab initio calcula- 
tions. In order to find the dynamical response within TD DFT, only the knowledge of these 
unperturbed Kohn-Sham states (both occupied and unoccupied) is required; no knowledge 
of real excitation spectra (both energies and lifetimes) is necessary. This is the main advan- 
tage of such approach. Unfortunately, within TD DFT, an accurate approximation to the 
kernel /xc(r, r', uj) describing dynamical exchange-correlation effects is unknown while some 
progress is currently been made []T9[. In the following, the static local density approximation 
(LDA) 1^ improved by a generalized gradient approximation (GGA) is adopted to treat 
/a;c(r, r', tu). To date, these are the most popular tools for practical ab initio calculations, 
which are known to produce static response functions as well as other ground-state, optical 
pl| , plus, recently [|lT|, superconducting and transport properties for large variety of solids 
in good agreement with experiments. The use of other approximations to J^c(r, r', tu) will 
be addressed in the future work. 

An important issue of variational linear-response formulation is now discussed. The 
advantage of variational principles for the calculation of physical quantities is that if one 
makes a first-order error in the trial function, the error in the variational quantity is of the 
second order. Static charge and spin susceptibilities appeared as second-order changes in 
the total energy due to applied external fields can be calculated in a variational way. This 
was demonstrated long time ago |^ on the example of magnetic response, and, recently 



pil| , |25[] , in the problem of lattice dynamics which is an example of charge response. The 
proof is directly related to a powerful "2n + 1" theorem of perturbation theory and station- 
arity property for the total energy itself |]2^. Any {2n + l)th change in the total energy 



Etot involves finding only (n)th order changes in one-electron wave functions ^pi, and corre- 
sponding changes in the charge density as well as in the magnetization. Any (2?T,)th change 
in Etot is then variational with respect to the (n)th-order changes in ipi. 

A time-dependent generalization of these results is now required. For TD external fields, 
the action S" as a functional of p{rt) and m(rt) is considered within TD DFT |T8| , p5| . These 



functions are expressed via Kohn-Sham spinor orbitals ipi (rt) satisfying TD Schrodinger's 
equation p^. Therefore, S as the stationary functional of ipi (r^) is considered in practice. 



When the external field is small, the perturbed wave function is represented as ipi (r) e -|- 
Sipi (rt) and the first-order changes (r^) define the induced charge density as well as the 
magnetization: 

6p = ^ {{6iJ,\I\fy + {Hl\5fy) (2) 

i 

Sm = pbY. {{^'^iWl'^i} + {^iWl^'^i}) (3) 

i 

Here { 1 1 } denotes averaging over spin degrees of freedom only, / is the unit 2x2 matrix, and 
a is the Pauli matrix. It is now seen that the knowledge of Sipt (rt) will solve the problem. 
In order to find Stpi (rt), a time-dependent analog of the "2n + 1" theorem is now in- 
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troduced. Any {2n + l)th change in the action functional 5* involves finding only (n)th 
order changes in the TD functions ipi (rt), and corresponding changes in charge density as 
well as in the magnetization. Any (2?T,)th change in S is then variational with respect to 
the (n)th-order changes in ipi (rt). The proof is the same as for the static case P3 if the 



stationarity property of 5* and the standard TD perturbation theory are exploited. For 
important case n = 2, this theorem makes the second-order change S**^^-* in the action vari- 
ational with respect to the first-order changes 6tpi (rt) . If the perturbation has the form 
(|^), S^'^^ is directly related to the real diagonal part of the dynamical spin susceptibility 
i?e[Xi/^(q + G', q + G,ci;)]g'=g ? thus allowing its variational estimate p5[] . 

The problem is now reduced to find S"*-^-* as a functional of 5ipi (vt) and to minimize 
it. This will bring an equation for 6tpi {rt). Any change in the action functional can be 
established by straightforward varying S of TD DFT |jl8,25| with respect to the perturbation 



(0). This is analogous to what is done in the static DFT to derive, for example, the dynamical 
matrix |T0|. S''-^^ is found to be 



i 

J 5p5V,ff - J 5m(5Be// + <5Be,,) (4) 

where the unperturbed 2x2 Hamiltonian matrix if = (— + V^ff)! — fiB<^^eff- Veff and 
Be// are the ground-state potential and magnetic field of the DFT. 5K// and ^Be// are 
their first-order changes induced by the perturbation (|lD which involve the Hartree (for 
SVeff) and the exchange-correlation contributions expressed via 5pand 6m in the standard 
manner 

The differential equation for Sipi (rt) is now derived from the stationarity condition of 
(I). It is given by 

(H - idtl)5^i + {SVeffI - fiB(r5B,ff)tf, = (5) 

This is a time-dependent version of the so-called Sternheimer equation which is the 
Schrodinger equation to linear order. It can be solved easily on the frequency axis which 
substitutes —idt by ei±uj in (H). The solution of the whole problem assumes self-consistency: 
First, Eq. is solved with the external field (|T]). Second, 6p{Tuj) and 5m(ra;) are found 
according to (0) and (0). Third, screened potential 6Veff{ru) and magnetic field 5Be//(rco') 
are constructed. The cycle is repeated again by solving (^. Evaluating S**-^-* after (^ brings 
the variational estimate of the real diagonal susceptibility at the iteration. The whole func- 
tion is accessed via the knowledge of 5m(ra;). The self-consistency should be done for every 
q + G and u value appeared in (|1|). 

The advantages of this method are now seen: First, Eq. (^ does not require an expansion 
of 6ipi over complete set of unperturbed wave functions ipj as it is done in the standard 
perturbation theory. Only the knowledge of occupied and those unoccupied states which 
are below Ep + a; is necessary. Second, the inversion problem is substituted by the self- 
consistent finding of 5K// and 5Be//. This normally requires about 10 iterations to reach the 
convergency. Third, the method treats on the same footing both longitudinal and transverse 
spin fiuctuations which is achieved by choosing the polarization 6h of the external filed (|l]) 
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along or perpendicular to the magnetization axe. Fourth, the method gives an access to 
c/iar^fg-spin fluctuations via the knowledge of 6p{ru!), and it is trivially converted to study 
dynamical charge fluctuations, if a TD scalar filed of the type is considered as the 
perturbation. 

An implementation of the method using linear muffin-tin orbital (LMTO) representation 
is now discussed. As the original wave function ipi is expanded in terms of the LMTOs Xa 
with the coefficients A^, the first-order change Sipi generally involves both changes 6A1 in 
the expansion coefficients and changes 6xa in the LMTO basis set ||10|- Changes 6 A"' are now 
new variational parameters instead of They must be found by minimizing the functional 
(I). Changes 6xa are, on the other hand, an auxiliary set of functions which is constructed 
to make the expansion of 6ipi fastly convergent. Basis {Sxa} is normally adjusted to the 
perturbation in the same way as the original basis {Xa} is tailored to the unperturbed one- 
electron potential. Such perturbative technique was found to be extremely efficient in the 
problem of lattice dynamics |jlO|. In the magnetic response calculation introducing 5xa is 
important for the fields exhibiting strong short-wavelength oscillations. On the other hand, 
in the calculations with G = in (P the contributions originating from 6xa are found to 
be small. 

Numerical efficiency of the method is now demonstrated by calculating spin suscepti- 
bilities for a number of metals. No shape approximations are made in these calculations 
either for the charge densities and the potentials or for the dynamical response functions. 
All the relevant quantities are expanded in spherical harmonics inside muffin-tin spheres 
and in plane waves in the interstitial region as it was done in original full-potential and 
static linear-response LMTO methods |TD[. The use of GGA for exchange and correlation 
gives practically coinciding theoretical and experimental lattice constants. Necessary Bril- 
louin zone (BZ) integrals are carried out using a multigrid tetrahedron technique with 
thousand k points. 

The ab initio results obtained for bcc Fe are now reported. Fig. |I] shows calculated 
transverse spin susceptibility /m[x+_(q, q, u)] for q = (00x)27r/a At small q the undecaying 
spin waves are seen to persist in the structure of exhibiting a standard dispersion law 

uj{q) = Dq^, where D is the spin stiffness of the material. The spin waves rapidly decay 
when q approaches to approximately one-half of the BZ. Similar picture has been found for 
the q's along (111) direction. The deduced spin-wave spectrum is shown by the solid line 
on top of Fig. 1. It agrees well with the experiment shown by circles as well as with the 
recent frozen-magnon calculations [|l^]. Also, in agreement with experiment any additional 
structure which can be attributed to the appearance of optical spin-wave branches is not 
predicted. This advances the early RPA calculation [0. 

Jm[x+_(q, q, cj)] obtained for fee Ni is shown on Fig. |^. The unusual structure for the 
energies near 100 meV and for the q's (0,0,0.2-0.4) 27r/a is clearly distinguishable. This 



was attributed to the appearance of the optical branch in the spin-wave spectrum |]7,I^ 
However, since this structure is seen to be only localized in a certain region of q space, its 
interpreting as a well-defined branch persisting to the BZ boundary might be not com- 
pletely correct. The computations along (111) direction do not show such unusual behavior. 
The obtained spin-wave spectrum (line on top of Fig. 2) is in agreement with the measured 
one (balls) [|1^] in the low-frequency interval. However, a tendency to overestimate spin- 
wave energies for higher u is found both for (001) and (111) directions. This is attributed 
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to the poor treatment of dynamical exchange-correlation effects due to simple GGA. 

Two examples of calculating paramagnetic spin fluctuations are now considered. Fig. |^ 
shows calculated Jm[x(q, q, u^)] for paramagnetic bcc Cr. A remarkable structure is clearly 
seen for the q's near {0,0,xsdw ~0.9)27r/a, where the susceptibility is mostly enhanced at 
low frequencies (experimentally, xsdw=0-Q5). This predicts Cr to be an incommensurate 
antiferromagnet. To clarify whether the Fermi-surface nesting is the origin of such behavior 
]T^ , the non-interacting susceptibility Im[xQ{(i,ci,uj)] can be analyzed. It does not show 
up a structure peaked at xsdw ~0.9, and is only a monotonically varying function when x 
increases from to 1. This means that the generalized Stoner criterium 1 = IxcXoi^) does 
not necessarily assumes a peak in Xoi^SDw) for Cr. 

/m[x(q, q, u)] in Pd is found to be strongly enhanced at small q's in complete agreement 
with the early studies |Q. Therefore, the method also confirms a closeness of Pd to the 
ferromagnetic instability. 

In conclusion, the developed approach is able to describe known spin-fiuctuational spec- 
tra of real materials which demonstrates its efficiency for practical ab initio calculations. 
Also, more elaborate approximations to the dynamical exchange and correlation are clearly 
required in order to account for the observed discrepancies. 

The author is indebted to O. K. Andersen, O. Gunnarsson, O. Jepsen, M. I. Katsnelson, 
and A. I. Liechtenstein for many helpful discussions. 
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FIGURES 



FIG. 1. Calculated Im[x-\ (q, q,a;)] (arb. units) for Fe. Top line shows the deduced magnon 

spectrum. Balls indicate the experimental data [14]. 

FIG. 2. Calculated Im[x-\ (q, q, w)] (arb. units) for Ni. Top line shows the deduced magnon 

spectrum. Balls indicate the experimental data [13]. 

FIG. 3. Calculated Jm[x(q,q,a;)] (Ry"^) for Cr. 
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Fig.l, Savrasov 




Fig. 2, Savrasov 




Fig. 3, Savrasov 




